#include<iostream>
#include<math.h>

double Getmax(double a,double b,double c)
{
    return std::max(a,std::max(b,c));
}

int main()
{
    int n=5;
    double f1[n+1][n+1],f2[n+1][n+1];
    int a;
    int count=0;//a 迭代判断 count 迭代次数
    double fai_n=200,fai_s=0,fai_w=0,fai_e=100,fai_i=0;//初值
    double rou=1.0;//密度
    double u=2.5;
    double v=2.5;//x方向速度 y方向速度
    double gama=0.1;//性质
    double len=1.0;
    double dx,dy;
    dx=dy=len/n;
    int i,j;
    double Fn,Fs,Fw,Fe;
    Fn=rou*v*dx;Fs=rou*v*dx;Fw=rou*u*dy;Fe=rou*u*dy;
    double Dn,Ds,Dw,De;
    Dn=gama*dx/dy;Ds=gama*dx/dy;Dw=gama*dy/dx;De=gama*dy/dx;

    for(i=0;i<=n+1;i++)//s
    {
	j=0;
	f1[i][j]=fai_s;
    }
    for(j=0;j<=n+1;j++)//w
    {
	i=0;
	f1[i][j]=fai_w;
    }
    for(i=0;i<=n+1;i++)//n
    {
	j=n+1;
	f1[i][j]=fai_n;
    }
    for(j=0;j<=n+1;j++)//e
    {
	i=n+1;
	f1[i][j]=fai_e;
    }
    for(i=1;i<n+1;i++)//内部
    {
	for(j=1;j<n+1;j++)
	{
	    f1[i][j]=fai_i;
	}
    }

    do
    {
	a=0;
	for(i=0;i<=n+1;i++)
	{
	    for(j=0;j<=n+1;j++)
	    {
		f2[i][j]=f1[i][j];
	    }
	}

	
	for(i=2;i<=n-1;i++)//内部
	{
	    for(j=2;j<=n-1;j++)
	    {
	        f1[i][j]=(Getmax(-Fn,(Dn-Fn/2),0)*f1[i][j+1]+Getmax(Fw,(Dw+Fw/2),0)*f1[i-1][j]+Getmax(Fs,(Ds+Fs/2),0)*f1[i][j-1]+Getmax(-Fe,(De-Fe/2),0)*f1[i+1][j])/(Getmax(Fw,(Dw+Fw/2),0)+Getmax(-Fe,(De-Fe/2),0)+Getmax(Fs,(Ds+Fs/2),0)+Getmax(-Fn,(Dn-Fn/2),0)+Fe-Fw+Fn-Fs);
	    }
	}
	for(i=2;i<=n-1;i++)//s
	{
	    j=1;
	    f1[i][j]=(Getmax(-Fn,(Dn-Fn/2),0)*f1[i][j+1]+Getmax(Fw,(Dw+Fw/2),0)*f1[i-1][j]+Getmax(Fs,(2*Ds+Fs/2),0)*f1[i][j-1]+Getmax(-Fe,(De-Fe/2),0)*f1[i+1][j])/(Getmax(Fw,(Dw+Fw/2),0)+Getmax(-Fe,(De-Fe/2),0)+Getmax(Fs,(2*Ds+Fs/2),0)+Getmax(-Fn,(Dn-Fn/2),0)+Fe-Fw+Fn-Fs);
	}
	for(j=2;j<=n-1;j++)//w
	{
	    i=1;
	    f1[i][j]=(Getmax(-Fn,(Dn-Fn/2),0)*f1[i][j+1]+Getmax(Fw,(2*Dw+Fw/2),0)*f1[i-1][j]+Getmax(Fs,(Ds+Fs/2),0)*f1[i][j-1]+Getmax(-Fe,(De-Fe/2),0)*f1[i+1][j])/(Getmax(Fw,(2*Dw+Fw/2),0)+Getmax(-Fe,(De-Fe/2),0)+Getmax(Fs,(Ds+Fs/2),0)+Getmax(-Fn,(Dn-Fn/2),0)+Fe-Fw+Fn-Fs);
	}
	for(i=2;i<=n-1;i++)//n
	{
	    j=n;
	    f1[i][j]=(Getmax(-Fn,(2*Dn-Fn/2),0)*f1[i][j+1]+Getmax(Fw,(Dw+Fw/2),0)*f1[i-1][j]+Getmax(Fs,(Ds+Fs/2),0)*f1[i][j-1]+Getmax(-Fe,(De-Fe/2),0)*f1[i+1][j])/(Getmax(Fw,(Dw+Fw/2),0)+Getmax(-Fe,(De-Fe/2),0)+Getmax(Fs,(Ds+Fs/2),0)+Getmax(-Fn,(2*Dn-Fn/2),0)+Fe-Fw+Fn-Fs);
	}
	for(j=2;j<=n-1;j++)//e
	{
	    i=n;
	    f1[i][j]=(Getmax(-Fn,(Dn-Fn/2),0)*f1[i][j+1]+Getmax(Fw,(Dw+Fw/2),0)*f1[i-1][j]+Getmax(Fs,(Ds+Fs/2),0)*f1[i][j-1]+Getmax(-Fe,(2*De-Fe/2),0)*f1[i+1][j])/(Getmax(Fw,(Dw+Fw/2),0)+Getmax(-Fe,(2*De-Fe/2),0)+Getmax(Fs,(Ds+Fs/2),0)+Getmax(-Fn,(Dn-Fn/2),0)+Fe-Fw+Fn-Fs);
	}
	i=1;j=1;//左下
	{
	    f1[i][j]=(Getmax(-Fn,(Dn-Fn/2),0)*f1[i][j+1]+Getmax(Fw,(2*Dw+Fw/2),0)*f1[i-1][j]+Getmax(Fs,(2*Ds+Fs/2),0)*f1[i][j-1]+Getmax(-Fe,(De-Fe/2),0)*f1[i+1][j])/(Getmax(Fw,(2*Dw+Fw/2),0)+Getmax(-Fe,(De-Fe/2),0)+Getmax(Fs,(2*Ds+Fs/2),0)+Getmax(-Fn,(Dn-Fn/2),0)+Fe-Fw+Fn-Fs);
	}
	i=1;j=n;//左上
	{
	    f1[i][j]=(Getmax(-Fn,(2*Dn-Fn/2),0)*f1[i][j+1]+Getmax(Fw,(2*Dw+Fw/2),0)*f1[i-1][j]+Getmax(Fs,(Ds+Fs/2),0)*f1[i][j-1]+Getmax(-Fe,(De-Fe/2),0)*f1[i+1][j])/(Getmax(Fw,(2*Dw+Fw/2),0)+Getmax(-Fe,(De-Fe/2),0)+Getmax(Fs,(Ds+Fs/2),0)+Getmax(-Fn,(2*Dn-Fn/2),0)+Fe-Fw+Fn-Fs);
	}
	i=n;j=1;//右下
	{
	    f1[i][j]=(Getmax(-Fn,(Dn-Fn/2),0)*f1[i][j+1]+Getmax(Fw,(Dw+Fw/2),0)*f1[i-1][j]+Getmax(Fs,(2*Ds+Fs/2),0)*f1[i][j-1]+Getmax(-Fe,(2*De-Fe/2),0)*f1[i+1][j])/(Getmax(Fw,(Dw+Fw/2),0)+Getmax(-Fe,(2*De-Fe/2),0)+Getmax(Fs,(2*Ds+Fs/2),0)+Getmax(-Fn,(Dn-Fn/2),0)+Fe-Fw+Fn-Fs);
	}
	i=n;j=n;//右上
	{
	    f1[i][j]=(Getmax(-Fn,(2*Dn-Fn/2),0)*f1[i][j+1]+Getmax(Fw,(Dw+Fw/2),0)*f1[i-1][j]+Getmax(Fs,(Ds+Fs/2),0)*f1[i][j-1]+Getmax(-Fe,(2*De-Fe/2),0)*f1[i+1][j])/(Getmax(Fw,(Dw+Fw/2),0)+Getmax(-Fe,(2*De-Fe/2),0)+Getmax(Fs,(Ds+Fs/2),0)+Getmax(-Fn,(2*Dn-Fn/2),0)+Fe-Fw+Fn-Fs);
	}


	for(i=1;i<n+1;i++)
	{
	    for(j=1;j<n+1;j++)
	    {
		if(fabs(f1[i][j]-f2[i][j])>=1e-3)
		{
		    a=1;
		}
	    }
	}
	count++;
    }while(a==1);

    for(j=1;j<n+1;j++)
    {
	for(i=1;i<n+1;i++)
	{
	    std::cout<<f1[i][j]<<" ";
	}
	std::cout<<std::endl;
    }
    std::cout<<"迭代次数="<<count<<std::endl;
}

